# ============================================================================
# Study 1 Replication Script: 
# Gender Differences in Immigration Attitudes in Japan 
# Starting from study1_replication_data.RData
# ============================================================================

# Load necessary libraries
library(tidyverse)  
library(broom)      
library(ggrepel)    
library(ggeffects)  
library(scales)     
library(stargazer)  
library(forcats)   

# ============================================================================
# LOAD REPLICATION DATA
# ============================================================================

# Load the replication dataset
load("study1_replication_data.RData")

# ============================================================================
# HELPER FUNCTIONS
# ============================================================================

# Function to calculate summary statistics with confidence intervals
summarySE <- function(data=NULL, measurevar, groupvars=NULL, na.rm=FALSE,
                      conf.interval=.95, .drop=TRUE) {
  
  length2 <- function (x, na.rm=FALSE) {
    if (na.rm) sum(!is.na(x))
    else       length(x)
  }
  
  datac <- plyr::ddply(data, groupvars, .drop=.drop,
                       .fun = function(xx, col) {
                         c(N    = length2(xx[[col]], na.rm=na.rm),
                           mean = mean   (xx[[col]], na.rm=na.rm),
                           sd   = sd     (xx[[col]], na.rm=na.rm)
                         )
                       },
                       measurevar
  )
  
  datac <- plyr::rename(datac, c("mean" = measurevar))
  datac$se <- datac$sd / sqrt(datac$N)
  ciMult <- qt(conf.interval/2 + .5, datac$N-1)
  datac$ci <- datac$se * ciMult
  
  return(datac)
}


# ============================================================================
# FIGURE 3: MEAN SCORES BY EXPERIMENTAL CONDITION AND GENDER (MANUSCRIPT)
# Note: This replication code produces a combined plot showing both outcomes
# (feelings toward immigrant women and men) with different line types.
# In the published manuscript, these are presented as separate Figure 3a and 3b
# for clearer visualization of each outcome.
# ============================================================================

# Prepare data for plotting
us_gender <- jp_replication[jp_replication$respondent_gender %in% c("Japanese Women", "Japanese Men"), ]

# Create combined data for both outcomes
combined_data <- bind_rows(
  us_gender %>%
    dplyr::group_by(conditions) %>%
    pivot_longer(cols = c('ft_wom_r'), names_to = "ft_var", values_to = "ft_score") %>%
    summarySE(measurevar = "ft_score", groupvars = c("respondent_gender", "conditions")) %>%
    mutate(outcome = "Feelings toward Immigrant Women"),
  
  us_gender %>%
    dplyr::group_by(conditions) %>%
    pivot_longer(cols = c('ft_men_r'), names_to = "ft_var", values_to = "ft_score") %>%
    summarySE(measurevar = "ft_score", groupvars = c("respondent_gender", "conditions")) %>%
    mutate(outcome = "Feelings toward Immigrant Men")
)

combined_data$conditions <- factor(combined_data$conditions, 
                                   levels = c("Control", "Treat1. Professional Immigration", "Treat2. Non-Professional Immigration"))
combined_data$outcome <- factor(combined_data$outcome, 
                                levels = c("Feelings toward Immigrant Women", "Feelings toward Immigrant Men"))

# Create Figure 3
figure3 <- ggplot(combined_data, aes(x = conditions, y = ft_score, fill = respondent_gender, 
                                     group = interaction(respondent_gender, outcome))) +
  geom_linerange(aes(ymin = ft_score - ci, ymax = ft_score + ci, color = respondent_gender, linetype = outcome), 
                 position = position_dodge(width = 1), size = 1.1) +
  geom_point(aes(shape = respondent_gender, color = respondent_gender), 
             position = position_dodge(width = 1), size = 6, alpha = 0.9) +
  facet_wrap(~ conditions, scales = "free_x", nrow = 1) +
  theme_bw() +
  theme(
    strip.text = element_text(size = 17, face = "bold"),
    strip.placement = "outside",  
    axis.text.y = element_text(size = 17, face = "bold"), 
    axis.title.y = element_text(size = 17, face = "bold"), 
    axis.text.x = element_blank(),
    axis.title.x = element_blank(),
    axis.ticks.x = element_blank(),
    legend.text = element_text(size = 17, face = "bold"),
    legend.position = "bottom",
    legend.box = "vertical"
  ) + 
  scale_shape_manual(values = c("Japanese Women" = 15, "Japanese Men" = 17)) + 
  scale_color_manual(values = c("Japanese Women" = "maroon", "Japanese Men" = "darkorange")) +
  scale_fill_manual(values = c("Japanese Women" = "maroon", "Japanese Men" = "darkorange")) +
  labs(y = "Feelings towards Immigrants", x = "Experimental Conditions") +
  guides(
    color = guide_legend(title = "Participant Gender"),
    shape = guide_legend(title = "Participant Gender"),
    linetype = guide_legend(title = "Outcome"),
    fill = "none"
  ) +
  scale_y_continuous(labels = scales::number_format(accuracy = 0.01))

print(figure3)

# ============================================================================
# FIGURE 4: Marginal Effect of Japanese Women’s Job Type on Feelings
# Note: In the published manuscript, these are presented as separate 
# Figure 4a and 4b for clearer visualization of each outcome.
# ============================================================================

# Run Regression Models
# Model 1: Feeling thermometer toward immigrant women
model_women <- lm(ft_wom_r ~ conditions:job.nonprof.binary:female +
                    conditions + job.nonprof.binary + female +
                    inc_11_r + edu_8_r + age_6_r + pride_national_r +
                    econ_nation_notsat_r + econ_self_notsat_r +
                    job_conf_na_mean + pid_ldp,
                  data = jp_replication)
# Model 2: Feeling thermometer toward immigrant men
model_men <- lm(ft_men_r ~ conditions:job.nonprof.binary:female +
                  conditions + job.nonprof.binary + female +
                  inc_11_r + edu_8_r + age_6_r + pride_national_r +
                  econ_nation_notsat_r + econ_self_notsat_r +
                  job_conf_na_mean + pid_ldp,
                data = jp_replication)

# Generate predicted values for immigrant women outcome
predictions_women <-  ggpredict(model_women, terms = c("conditions", "job.nonprof.binary", "female")) %>%
  filter(facet == 1) %>%  # Japanese women only
  mutate(outcome = "Feelings toward Immigrant Women")

# Generate predicted values for immigrant men outcome  
predictions_men <-  ggpredict(model_men, terms = c("conditions", "job.nonprof.binary", "female")) %>%
  filter(facet == 1) %>%  # Japanese women only
  mutate(outcome = "Feelings toward Immigrants Men")

# Combine predictions from both models for comparative visualization
marginal_effects <- bind_rows(predictions_women, predictions_men) %>%
  filter(!is.na(x)) # Extra step to filter out NAs created by ggpredict()

# Set factor levels for proper ordering in plot
marginal_effects$x <- factor(marginal_effects$x, levels = c("Control", 
                                                            "Treat1. Professional Immigration", 
                                                            "Treat2. Non-Professional Immigration"))

# Set outcome factor levels to match the correct legend order
marginal_effects$outcome <- factor(marginal_effects$outcome, 
                                   levels = c("Feelings toward Immigrant Women", "Feelings toward Immigrants Men"))

marginal_effects <- marginal_effects %>%
  filter(!is.na(x)) # Extra step to filter out NAs created by ggpredict()

# Create Figure 4: Marginal effects plot with predicted values and confidence intervals
figure4 <- ggplot(marginal_effects, aes(y = predicted, x = factor(x), 
                                        color = factor(group), shape = factor(group), 
                                        linetype = outcome, group = interaction(group, outcome))) +
  geom_point(position = position_dodge(width = 1), size = 5, na.rm = TRUE) +
  geom_linerange(aes(ymin = conf.low, ymax = conf.high), 
                 position = position_dodge(width = 1), linewidth = 1.1, na.rm = TRUE) +
  labs(x = "Participant's Job Type",
       y = "Predicted Values of Japanese Women's Feelings",
       color = "Participant's Job Type",
       shape = "Participant's Job Type",
       linetype = "Outcome") + 
  scale_color_manual(values = c("forestgreen", "darkslateblue"), 
                     labels = c("Professional Jobs", "Non-Professional Jobs")) +
  scale_shape_manual(values = c(16, 16), 
                     labels = c("Professional Jobs", "Non-Professional Jobs")) +
  scale_linetype_manual(values = c("solid", "dotted")) +
  facet_wrap(~ x, scales = "free_x", nrow = 1) +
  theme_bw() +
  theme(strip.text = element_text(size = 17, face = "bold"),
        axis.title.y = element_text(size = 17, face = "bold"),
        axis.text.x = element_blank(),    
        axis.title.x = element_blank(),   
        axis.ticks.x = element_blank(),   
        axis.text.y = element_text(size = 17, face = "bold"),
        legend.text = element_text(size = 17, face = "bold"),
        plot.title = element_text(size = 20, face = "bold", hjust = 0.5),
        legend.position = "bottom",
        legend.box = "vertical") + 
  guides(color = guide_legend(title = "Participant's Job Type"),
         shape = guide_legend(title = "Participant's Job Type"),
         linetype = guide_legend(title = "Outcome")) +
  scale_y_continuous(labels = scales::number_format(accuracy = 0.01))

print(figure4)

# ============================================================================
# FIGURE 5: Marginal Effect of Japanese Men’s Job Type on Feelings
# Note: In the published manuscript, these are presented as separate 
# Figure 5a and 5b for clearer visualization of each outcome.
# ============================================================================

# Run Regression Models
# Model 1: Feeling thermometer toward immigrant women
model_women <- lm(ft_wom_r ~ conditions:job.nonprof.binary:female +
                    conditions + job.nonprof.binary + female +
                    inc_11_r + edu_8_r + age_6_r + pride_national_r +
                    econ_nation_notsat_r + econ_self_notsat_r +
                    job_conf_na_mean + pid_ldp,
                  data = jp_replication)
# Model 2: Feeling thermometer toward immigrant men
model_men <- lm(ft_men_r ~ conditions:job.nonprof.binary:female +
                  conditions + job.nonprof.binary + female +
                  inc_11_r + edu_8_r + age_6_r + pride_national_r +
                  econ_nation_notsat_r + econ_self_notsat_r +
                  job_conf_na_mean + pid_ldp,
                data = jp_replication)

# Generate predicted values for immigrant women outcome
predictions_women <- ggpredict(model_women, terms = c("conditions", "job.nonprof.binary", "female")) %>%
  filter(facet == 0) %>%  # Japanese men only
  mutate(outcome = "Feelings toward Immigrant Women") %>%
  filter(!is.na(x)) # Extra step to filter out NAs created by ggpredict()

# Generate predicted values for immigrant men outcome  
predictions_men <- ggpredict(model_men, terms = c("conditions", "job.nonprof.binary", "female")) %>%
  filter(facet == 0) %>%  # Japanese men only
  mutate(outcome = "Feelings toward Immigrant Men") %>%
  filter(!is.na(x)) # Extra step to filter out NAs created by ggpredict()

# Combine predictions from both models for comparative visualization
marginal_effects_men <- bind_rows(predictions_women, predictions_men) %>%
  filter(!is.na(x)) # Extra step to filter out NAs created by ggpredict()

# Set factor levels for proper ordering in plot
marginal_effects_men$x <- factor(marginal_effects_men$x, levels = c("Control", 
                                                            "Treat1. Professional Immigration", 
                                                            "Treat2. Non-Professional Immigration"))

# Set outcome factor levels to match the correct legend order
marginal_effects_men$outcome <- factor(marginal_effects_men$outcome, 
                                   levels = c("Feelings toward Immigrant Women", "Feelings toward Immigrant Men"))

marginal_effects_men <- marginal_effects_men %>%
  filter(!is.na(x)) # Extra step to filter out NAs created by ggpredict()

# Create Figure 5: Marginal effects plot for Japanese men with predicted values and confidence intervals
figure5 <- ggplot(marginal_effects_men, aes(y = predicted, x = factor(x), 
                                            color = factor(group), shape = factor(group), 
                                            linetype = outcome, group = interaction(group, outcome))) +
  geom_point(position = position_dodge(width = 1), size = 5, na.rm = TRUE) +
  geom_linerange(aes(ymin = conf.low, ymax = conf.high), 
                 position = position_dodge(width = 1), linewidth = 1.1, na.rm = TRUE) +
  labs(x = "Participant's Job Type",
       y = "Predicted Values of Japanese Men's Feelings",
       color = "Participant's Job Type",
       shape = "Participant's Job Type",
       linetype = "Outcome") + 
  scale_color_manual(values = c("#6B8E23", "darkblue"), 
                     labels = c("Professional Jobs", "Non-Professional Jobs")) +
  scale_shape_manual(values = c(16, 16), 
                     labels = c("Professional Jobs", "Non-Professional Jobs")) +
  scale_linetype_manual(values = c("solid", "dotted")) +
  facet_wrap(~ x, scales = "free_x", nrow = 1) +
  theme_bw() +
  theme(strip.text = element_text(size = 17, face = "bold"),
        axis.title.y = element_text(size = 17, face = "bold"),
        axis.text.x = element_blank(),    # Hide x-axis labels for cleaner appearance
        axis.title.x = element_blank(),   # Condition names shown in facet strips
        axis.ticks.x = element_blank(),   # Remove x-axis tick marks
        axis.text.y = element_text(size = 17, face = "bold"),
        legend.text = element_text(size = 17, face = "bold"),
        plot.title = element_text(size = 20, face = "bold", hjust = 0.5),
        legend.position = "bottom",
        legend.box = "vertical") + 
  guides(color = guide_legend(title = "Participant's Job Type"),
         shape = guide_legend(title = "Participant's Job Type"),
         linetype = guide_legend(title = "Outcome")) +
  scale_y_continuous(labels = scales::number_format(accuracy = 0.01))

print(figure5)


# ============================================================================
# APPENDIX C2: JOB DISTRIBUTION BY GENDER
# ============================================================================
library(dplyr)
library(tidyr)
library(ggplot2)
library(forcats)

job_counts <- jp_replication %>%
  count(job, female) %>%
  pivot_wider(
    names_from = female,
    values_from = n,
    values_fill = list(n = 0),
    names_prefix = "gender_"
  ) %>%
  rename(male = gender_0, female = gender_1) %>%
  mutate(total = male + female) %>%
  arrange(desc(total)) %>%
  mutate(
    job_category = case_when(
      job == "1. 無職" ~ "Unemployed",
      job == "2. 学生" ~ "Student",
      job == "3. 専業主婦" ~ "Homemaker",
      job == "4. 退職した" ~ "Retired",
      job == "5. 自営業" ~ "Self-employed",
      job == "6. 事務職" ~ "Clerical",
      job == "7. 経営/管理職" ~ "Management/administrative",
      job == "8. 専門職/準専門職" ~ "(Semi-) Professional",
      job == "9. 熟練技能職" ~ "(Semi-) Skilled",
      job == "10. 販売/サービス業" ~ "Sales/Service",
      job == "11. 農林畜水産業" ~ "Agriculture/forestry/livestock/fishery",
      job == "12. 生産/運輸/日用/作業職" ~ "Production/Transportation/Labor",
      job == "13. 公務員" ~ "Government employee",
      job == "14. 軍人" ~ "Military",
      job == "15. その他 (答えを記入してください)" ~ "Other",
      TRUE ~ as.character(job)
    ),
    job_category = fct_reorder(job_category, total)  # reorder by total
  )

# Reshape for stacked bar chart
job_data_long <- job_counts %>%
  pivot_longer(cols = c(male, female), names_to = "gender", values_to = "count") %>%
  dplyr::select(job_category, gender, count)

# Plot
figure_c2 <- ggplot(job_data_long, aes(x = count, y = job_category, fill = gender)) +
  geom_bar(stat = "identity", position = "stack") +
  geom_text(aes(label = ifelse(count > 0, count, "")),
            position = position_stack(vjust = 0.5),
            color = "black", fontface = "bold", size = 3.5) +
  scale_fill_manual(values = c("male" = "#4A90E2", "female" = "#E24A90"),
                    labels = c("Male", "Female")) +
  labs(x = "Count", y = "Job Categories", fill = "Gender") +
  theme_minimal() +
  theme(axis.text.y = element_text(size = 11, face="bold"),
        legend.position = "bottom",
        legend.text = element_text(size = 12)) +
  scale_x_continuous(breaks = seq(0, 250, 50))

print(figure_c2)

# ============================================================================
# APPENDIX FIGURES D1 & D2: POST-TREATMENT ANALYSIS
# ============================================================================

# Appendix D1: Gender of Immigrant Thought of
prop_prof <- 100 * prop.table(table(jp_replication$imagine_sex_h))
prop_nonp <- 100 * prop.table(table(jp_replication$imagine_sex_l))

df_img <- data.frame(
  category = c("Female", "Male", "Other", "Female", "Male", "Other"),
  proportion = c(prop_prof["Female"], prop_prof["Male"], prop_prof["Other"], 
                 prop_nonp["Female"], prop_nonp["Male"], prop_nonp["Other"]),
  dataset = c(rep("Treat1.Professional Immigrants", 3), rep("Treat2.Non-Professional Immigrants", 3))
)

df_img$category <- factor(df_img$category, levels = c("Male", "Female", "Other"))
df_img$proportion_label <- sprintf("%.1f%%", df_img$proportion)

# Updated colors to match Image 1: Orange for Male, Teal for Female, Purple for Other
fill_colors <- c("#FF7F50", "#4ECDC4", "#9B9AE8")  # Orange, Teal, Purple

# Appendix Figure D1: Gender of Immigrant Thought of
figure_d1 <- ggplot(df_img, aes(x = dataset, y = proportion, fill = category)) +
  geom_bar(stat = "identity", position = "dodge", width = 0.7) +
  geom_text(aes(label = proportion_label),
            position = position_dodge(width = 0.7), vjust = -0.5, size = 3.5) +
  labs(
    x = "Experimental Conditions",
    y = "Percentage",
    fill = "Immigrant's Gender Thought"
  ) +
  scale_fill_manual(values = fill_colors) +
  theme_minimal() +
  theme(
    axis.title.y = element_text(size = 13),
    axis.text.x = element_text(size = 11, face = "bold"), 
    axis.title.x = element_text(size = 13), 
    legend.text = element_text(size = 13)
  )

print(figure_d1)

# Appendix D2: Origin of Immigrant Thought of
prop_table_t1 <- 100 * prop.table(table(jp_replication$imagine.ori.h))
prop_table_t2 <- 100 * prop.table(table(jp_replication$imagine.ori.l))

df_t1 <- data.frame(
  origins_thought = names(prop_table_t1),
  proportion = as.numeric(prop_table_t1),
  dataset = "Treat1.Professional Immigrants"
)

df_t2 <- data.frame(
  origins_thought = names(prop_table_t2),
  proportion = as.numeric(prop_table_t2),
  dataset = "Treat2.Non-Professional Immigrants"
)

df_combined <- rbind(df_t1, df_t2)

# Add "Eurasia" with 0% proportion if missing
if(!"Eurasia" %in% df_combined$origins_thought[df_combined$dataset == "Treat2.Non-Professional Immigrants"]) {
  df_combined <- rbind(df_combined,
                       data.frame(origins_thought = "Eurasia",
                                  proportion = 0,
                                  dataset = "Treat2.Non-Professional Immigrants"))
}

combined_proportions <- aggregate(proportion ~ origins_thought, data = df_combined, sum)
df_combined$origins_thought <- factor(df_combined$origins_thought, 
                                      levels = combined_proportions$origins_thought[order(combined_proportions$proportion, decreasing = TRUE)])

df_combined$proportion_label <- sprintf("%.1f%%", df_combined$proportion)

# Updated colors to match Images 2 & 3: Vibrant colors for different origins
fill_colors_origin <- c(
  "#4ECDC4",  # Teal for Southeast Asia (highest)
  "#FF7F50",  # Orange for East Asia  
  "#F39C12",  # Yellow for South Asia
  "#D2B48C",  # Tan for North America
  "#6495ED",  # Blue for Europe
  "#808080",  # Gray for Africa
  "#B19CD9",  # Light Purple for Other
  "#9ACD32",  # Green for Middle East
  "#FF69B4",  # Pink for Latin America
  "#FF8C00"   # Dark Orange for Eurasia
)

# Appendix Figure D2: Origin of Immigrant Thought of
figure_d2 <- ggplot(df_combined, aes(x = proportion, y = fct_reorder(origins_thought, proportion), fill = origins_thought)) +
  geom_bar(stat = "identity", position = position_dodge(width = 0.9), width = 0.4) +
  geom_text(aes(label = proportion_label),
            position = position_dodge(width = 0.9), hjust = -0.1, 
            size = 3, fontface = "bold") +
  facet_wrap(~dataset, scales = "free_y") +
  labs(
    x = "Percentage",
    y = "Immigrant's Origin Thought"
  ) +
  scale_fill_manual(values = fill_colors_origin) +
  theme_minimal() +
  theme(
    axis.title.y = element_text(size = 13),
    axis.text.y = element_text(size = 13), 
    axis.text.x = element_text(size = 13), 
    axis.title.x = element_text(size = 13),
    legend.position = "none",
    strip.text = element_text(size = 13, face = "bold"),
    plot.margin = margin(t = 1, r = 0.5, b = 1, l = 0.5, unit = "cm")
  )

print(figure_d2)


# ============================================================================
# APPENDIX E: REGRESSION RESULTS TABLES (F1 & F2)
# ============================================================================

library(stargazer)

# Appendix E1: Results for Feelings toward Immigrant Women (OLS)
model_women_base <- lm(ft_wom_r ~ conditions + job.nonprof.binary + female +
                         inc_11_r + edu_8_r + age_6_r + pride_national_r +
                         econ_nation_notsat_r + econ_self_notsat_r +
                         job_conf_na_mean + pid_ldp,
                       data = jp_replication)

stargazer(model_women_base, model_women, 
          title="",
          header=FALSE, 
          single.row = TRUE, 
          no.space = TRUE, 
          font.size = "small",
          dep.var.labels = "Feelings toward Immigrant Women",
          order = c("conditionsTreat2. Non-Professional Immigration:job.nonprof.binary:female",
                    "conditionsTreat1. Professional Immigration:job.nonprof.binary:female", 
                    "conditionsControl:job.nonprof.binary:female",
                    "conditionsTreat2. Non-Professional Immigration",
                    "conditionsTreat1. Professional Immigration", 
                    "job.nonprof.binary",
                    "female",
                    "inc_11_r",
                    "edu_8_r", 
                    "age_6_r",
                    "pride_national_r",
                    "econ_nation_notsat_r",
                    "econ_self_notsat_r",
                    "job_conf_na_mean",
                    "pid_ldp"),
          covariate.labels = c("Treat2:Job(NonProf):female",
                               "Treat1:Job(NonProf):female",
                               "Control:Job(NonProf):female", 
                               "Treat2(NonProf.Immi)",
                               "Treat1(Prof.Immi)",
                               "Participants.Job(NonProf)",
                               "Female", 
                               "Income",
                               "Education",
                               "Age",
                               "Nationality.Pride",
                               "National.Economy.Dissatisfied", 
                               "Own.Finance.Dissatisfied",
                               "Job.Security",
                               "pid.LDP"),
          flip = TRUE)

# Appendix E2: Results for Feelings toward Immigrant Men (OLS)
model_men_base <- lm(ft_men_r ~ conditions + job.nonprof.binary + female +
                         inc_11_r + edu_8_r + age_6_r + pride_national_r +
                         econ_nation_notsat_r + econ_self_notsat_r +
                         job_conf_na_mean + pid_ldp,
                       data = jp_replication)

stargazer(model_men_base, model_men, 
          title="",
          header=FALSE, 
          single.row = TRUE, 
          no.space = TRUE, 
          font.size = "small",
          dep.var.labels = "Feelings toward Immigrant Men",
          order = c("conditionsTreat2. Non-Professional Immigration:job.nonprof.binary:female",
                    "conditionsTreat1. Professional Immigration:job.nonprof.binary:female", 
                    "conditionsControl:job.nonprof.binary:female",
                    "conditionsTreat2. Non-Professional Immigration",
                    "conditionsTreat1. Professional Immigration", 
                    "job.nonprof.binary",
                    "female",
                    "inc_11_r",
                    "edu_8_r", 
                    "age_6_r",
                    "pride_national_r",
                    "econ_nation_notsat_r",
                    "econ_self_notsat_r",
                    "job_conf_na_mean",
                    "pid_ldp"),
          covariate.labels = c("Treat2:Job(NonProf):female",
                               "Treat1:Job(NonProf):female",
                               "Control:Job(NonProf):female", 
                               "Treat2(NonProf.Immi)",
                               "Treat1(Prof.Immi)",
                               "Participants.Job(NonProf)",
                               "Female", 
                               "Income",
                               "Education",
                               "Age",
                               "Nationality.Pride",
                               "National.Economy.Dissatisfied", 
                               "Own.Finance.Dissatisfied",
                               "Job.Security",
                               "pid.LDP"),
          flip = TRUE)

